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ABSTRACT 

A  distinction  is  often  in;ule  between  domain  decomposition  algorillims  for  elliptic  partial 
dilTcrential  equations  which  use  overlapping  subrcgions  and  those  which  do  not. 
Schwar/.'s  alternating  methoil,  the  oldest  of  them  all,  belongs  to  the  first  category.  It  has 
recently  been  discovered  that  many  of  the  algorithms  that  belong  to  the  second  category 
can  also  be  regarded  as  generalizations  of  the  classical  Schwar/  mettiod  or  nn  additive 
variant  thereof.  This  neu  approach  provides  new  tools  for  the  analysis  and  development 
of  domain  decomposition  algorithms.  In  this  paper,  we  introduce  an  abstract  additive 
Schwarz  method  and  de\elop  a  simple  framework  for  the  analysis  of  its  rate  of  conver- 
gence. We  show  that  it  is  pc^ssiblc  and  convenient  to  specify  algorithms  in  terms  of  a  set 
of  subspaces  and  rel;itcd  orthogonal  projections.  This  family  of  algorithms  is  further  ex- 
panded by  replacing  the  linear  systems  of  equations,  whicti  correspond  to  the  projec- 
tions, by  suitable  precoinlili<nuTS. 

We  illustrate  the  usefulness  of  this  approach  by  considering  additive  Schwar/  methods  of 
a  relatively  conventional  type,  iterative  substrucluring  methods,  domain  decomposition 
methods  defined  on  the  curves  or  surfaces  whicli  subdivide  the  region,  and  iterative 
refinement  methods.  Throughout  wc  work  in  a  framework  of  confoiming  finite  eleinents 
and  self-adjoint  problems,  but  we  also  mention  some  new  results  for  more  general  ellip- 
tic finite  element  pioblcms. 

Key  Words  Schw,ir/'s  altcrn.iting  method,  domain  dccom|->osilion,  elliptic  equations, 
finite  elements. 
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1.  Introduction.  In  this  paper,  we  survey  our  recent  research  on  domain  de- 
composition and  related  algorithms  for  elliptic  finite  element  problems.  One  of  our 
aims  is  to  develop  a  set  of  common  and  powerful  tools  for  the  development  and  anal- 
ysis of  a  great  variety  of  methods.  We  have  recently  discovered  that  many  iterative 
substructuring  methods  can  be  viewed  as  so-called  additive  Schwarz  methods  and 
that  the  analysis  of  these  methods  can  be  simplified;  cf.  Dryja  and  Widlund  [28]. 
Here,  we  discuss  and  extend  this  result  and  provide  further  evidence  that  many  do- 
main decomposition  methods  fit  into  a  framework  provided  by  the  Schwarz  additive 
methods. 

We  note  that  the  bounds  given  for  the  condition  number  of  the  iteration  operators 
for  many  domain  decomposition  algorithms  considered  in  this  and  other  papers  are 
either  uniform  in  the  number  of  subdomains  and  subregions  or  grow  only  polynomially 
in  the  logarithm  of  the  number  of  degrees  of  freedom  associated  with  an  individual 
subregion.  The  algorithms  are  therefore,  in  a  certain  sense,  almost  optimal. 

The  paper  is  organized  as  follows.  After  introducing  two  elliptic  model  problems 
and  certain  finite  element  methods  in  Section  2,  we  begin  Section  3  by  reviewing 
Schwarz's  alternating  algorithm  in  its  classical  setting.  Following  Sobolev  [50]  and  P. 
L.  Lions  [34],  we  indicate  how  this  algorithm  can  be  expressed  in  a  variational  form. 
Since  this  formulation  is  very  convenient  for  the  analysis  of  finite  element  problems, 
we  work  in  this  Hilbert  space  setting  throughout  the  paper. 

In  Section  3,  we  also  introduce  the  additive  variant  of  Schwarz's  method;  cf  Dryja 
[21],  Dryja  and  Widlund  [26],  [28],  Matsokin  and  Nepomnyaschikh  [41]  and  Nepom- 
nyaschikh  [44].  We  present  it  in  a  general  form  and  show  how  algorithms  of  this  kind 
can  be  defined  in  terms  of  a  set  of  subspaces  and  projections.  An  additive  Schwarz 
method  can  be  viewed  as  an  iterative  method  for  the  solution  of  an  auxiliary  linear 
problem  that  has  the  same  solution  as  the  original  finite  element  problem.  We  also 
introduce  tools,  which  make  it  possible  to  estimate  the  condition  number  k(P)  of  the 
operator  P  of  this  new  problem. 

It  is  often  straightforward  to  obtain  an  upper  bound  for  the  spectrum  of  P.  A 
lower  bound  of  the  eigenvalues  is  obtained  in  terms  of  an  upper  bound  of  a  Rayleigh 
quotient  which  measures  the  extent  by  which  the  subspaces  are  linearly  independent. 
(If  the  subspaces  are  orthogonal,  the  Schwarz  method  converges  in  one  step.)  We 
show  that  preconditioners  can  be  used  as  approximate  solvers  for  the  linear  systems, 
which  are  related  to  the  individual  projections,  and  that  it  is  easy  to  obtain  bounds 
for  the  condition  number  of  the  resulting  algorithm  in  terms  of  the  spectrum  of  the 
basic  operator  and  any  bounds  for  the  individual  preconditioned  problems  that  might 
be  available. 

Section  3,  thus,  provides  us  with  a  basic  framework.  In  the  rest  of  the  paper,  we 
turn  to  a  series  of  applications. 

In  Section  4,  we  summarize  our  work  on  additive  Schwarz  methods  of  a  type 
directly  related  to  the  classical  case.  For  previous  discussions  of  this  algorithm,  see 
Dryja  [21],  Dryja  and  Widlund  [26],  [28],  with  the  most  details  given  in  [28].  We  note 
that  this  method  has  been  extended  successfully  to  certain  stationary  and  parabolic 
convection-diffusion  problems  in  the  dissertation  of  Xiao-Chuan  Cai  [20], [19],  and  to 
mixed  finite  element  methods  of  the  Raviart-Thomas  type,  cf.  [46],  in  the  dissertation 
of  Tarek  Mathew  [40].  In  this  paper,  we  are  unable  to  provide  details  about  their 
work  since  it  would  require  the  introduction  of  many  notations,  etc.  We  also  discuss 
a  result  on  multilevel  algorithms  obtained  in  collaboration  with  Xuejun  Zhang,  a 
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Courant  Institute  graduate  student. 

In  Section  5,  we  discuss  a  well-known  family  of  domain  decomposition  methods 
known  as  iterative  substructuring  methods;  cf.  e.g.  Agoshkov  [1],  Bj0rstad  and  Hvid- 
sten  [6],  Bjorstad  and  Widlund  [8],  Bramble,  Pasciak  and  Schatz  [18],  [15],  [16],  [14], 
[17],  Dryja  [24],  Dryja,  Proskurowski  and  Widlund  [25],  Dryja  and  Widlund  [28],  Lebe- 
dev  [32],  Marchuk,  Kuznetsov  and  Matsokin  [39],  Quarteroni  [45],  Smith  and  Widlund 
[49]  and  Widlund  [54].  These  algorithms  are  ba.sed  on  a  non-overlapping  subdivision 
{fi,  }  of  the  region  fi.  (  Borrowing  a  term  from  structural  engineering,  the  subregions 
are  often  called  substructures.)  Surprisingly  enough,  as  shown  already  in  Dryja  and 
Widlund  [28],  many  of  these  methods  also  fit  well  in  the  Schwarz  framework.  In  this 
paper,  we  further  extend  our  analysis  and  include  a  more  detailed  discussion  of  the 
three-dimensional  case.  In  this  study,  it  is  fruitful  to  take  a  substructure  by  sub- 
structure view,  estimating  the  contribution  to  the  strain  energy  from  an  individual 
substructure  in  terms  of  the  corresponding  local  contribution  to  the  preconditioner. 
This  approach  is  similar  to  that  of  Jan  Mandel,  and  his  coworkers  Babuska,  Craig  and 
Pitkaranta,  cf.  [2],  [37],  [36],  [35],  who  have  begun  a  systematic  study  of  domain  de- 
composition methods  for  p-version  finite  element  methods;  cf.  also  Bramble,  Pasciak 
and  Schatz  [17]  for  a  fundamental  study  of  h-version  methods.  This  approach  high- 
lights how  preconditioners  can  be  constructed  from  parts  which  are  strictly  local  with 
respect  to  an  individual  substructure,  parts  which  involve  interaction  between  pairs 
of  neighboring  substructures  and  a  coarse  global  model  with  relatively  few  degrees  of 
freedom.  Many  of  the  results  which  were  obtained  for  two  subregions  in  the  early  de- 
velopment of  the  theory  can  now  be  recycled  to  yield  useful  results  for  the  much  more 
interesting  ca^e  of  many  substructures.  We  note  that  some  of  our  results  have  been 
extended  to  a  class  of  two-dimensional,  stationary  and  parabolic  convection-diffusion 
problems  by  Xiao-Chuan  Cai  [20],  [19]. 

In  Section  6,  we  consider  Schwarz  methods  on  the  lower  dimensional  manifolds 
formed  by  the  curves  or  surfaces,  which  partition  the  domain  into  substructures.  These 
algorithms  can  also  be  viewed  as  Sch%varz  methods  for  problems  of  potential  theory. 
This  idea  has  previously  been  discussed  by  Nepomnyaschikh  [44]  for  the  case  of  a  few 
substructures  and  without  a  full  development  of  the  theory.  We  include  a  discussion 
and  estimates  of  the  rate  of  convergence  of  a  method  developed  by  Bourgat,  Glowinski, 
Le  Tallec  and  Vidrascu  [11]  and  also  discuss  some  recent  work  by  Barry  Smith,  a 
Courant  Institute  graduate  student;  cf.  [48]. 

In  Section  7,  we  survey  our  recent  work  on  iterative  refinement  methods.  These 
are  methods  for  the  solution  of  the  linear  systems  of  algebraic  equations,  which  arise 
from  elliptic  finite  element  problems  defined  on  composite  meshes.  The  study  of  these 
methods  was  pioneered  by  McCormick  and  his  coworkers  Hart  and  Thomas  [31],  [42], 
[43].  An  analysis  of  two-level  algorithms  is  given  in  Bramble,  Ewing,  Pasciak  and 
Schatz  [13]  and  in  Mandel  and  McCormick  [38].  Proofs  of  the  results  discussed  in  this 
section  are  given  in  Dryja  and  Widlund  [27]  and  Widlund  [56].  The  main  difficulties  in 
this  work  are  related  to  an  effort  of  obtaining  bounds  for  the  rate  of  convergence  which 
are  independent  of  the  number  of  refinement  levels  as  well  as  the  number  of  degrees 
of  freedom.  We  note  that  in  the  dissertation  of  Tarek  Mathew  [40],  certain  results  are 
obtained  on  iterative  refinement  methods  for  Raviart-Thomas  finite  elements. 

In  this  paper,  we  primarily  consider  algorithms  that  are  of  additive  Schwarz  type. 
We  note  that  a  variety  of  multiplicative  algorithms  also  can  be  defined  systematically. 
At  present,  the  general  multiplicative  case  appears  to  be  less  well  understood  than 
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the  additive,  except  in  the  case  of  two  subspaces.  For  the  two  subspace  ca^e,  Bjorstad 
and  Mandel  [7]  have  recently  extended  their  earlier  work,  cf.  [5]  and  [38],  obtaining  a 
detailed  comparison  of  the  spectra  of  the  additive  and  multiplicative  algorithms.  We 
note  that  there  are  cases  for  which  a  satisfactory  theory  already  exists  in  the  case  of 
more  than  two  subspaces.  Thus  the  algorithm  of  Bank,  Dupont  and  Yserentant  [3]  is  a 
multiplicative  variant  of  the  original,  additive  algorithm  of  Yserentant's  [58];  it  is  just 
as  well  understood.  Similarly,  the  basic  theory  for  the  multilevel  iterative  refinement 
methods  has  advanced  to  the  same  level  in  the  two  cases;  cf.  Widlund  [56]  and  Dryja 
and  Widlund  [27].  We  note,  finally,  that  in  the  recent  thesis  of  Mathew  important 
progress  is  reported  for  the  general  multiplicative  case;  cf.  Mathew  [40]  and  a  remark 
in  Section  3. 

2.  Model  Problems  and  Finite  Elements.  In  this  section,  we  introduce  finite 
element  approximations  of  a  standard  Poisson  equation  and  a  special  second  order 
elliptic  problem  with  variable  coefficients. 

In  the  model  problems,  the  continuous  and  discrete  problems  are  of  the  form 

a{u,v)  =  f{v),    V  f  G    V, 

and 

(1)  a(uh,Vh)  =  fivh),    Vvh  €    V\ 

respectively.  We  consider  homogeneous  Dirichlet  problems  on  bounded  Lipschitz  re- 
gions in  two  or  three  dimensions  and  continuous,  piecewise  linear  finite  elements.  For 
the  Poisson  problem,  the  bilinear  form  is  defined  by 


(2)  a(u,v)=    f 


Vu  •  Vu  dx  . 

This  form  defines  a  semi-norm  |u|^wQx  =  (a(u,u))^/^  in  the  Sobolev  space  .ff^(f2).  It 
is  a  norm  of  V  =  HK^).  Here  Hq(Q)  is  the  subspace  of  H^{Q)  functions  with  zero 
trace;  aU  elements  of  V  and  its  subspace  V^  vanish  on  dfl,  the  boundary  of  fi.  We 
note  that  if  fi  C  0,  then  any  element  of  Hq{Ci)  can  be  extended  by  zero  to  an  element 
in  Hq{Q,),  and  that  the  extension  operator  is  continuous.  We  therefore  regard  Hq{CI) 
as  a  subspace  of  ^o(^)- 

Almost  all  our  results  can  be  extended  immediately  to  general  conforming  finite 
element  approximations  of  any  self-adjoint  elliptic  problem,  which  can  be  formulated 
as  a  minimization  problem,  and  to  more  general  boundary  conditions.  However,  some 
of  the  bounds  can  be  very  poor  if  there  is  a  great  variation  in  the  values  of  the 
coefficients.  In  this  paper,  we  consider  only  one  generalization, 

(3)  a{u,v)=   /  p{x)'V u  •  V V  dx  , 

Jn 

where  p(x)  >  0  can  be  discontinuous,  but  varies  slowly  in  each  subregion.  There  is 
experimental  evidence  that  some  domain  decomposition  methods  of  Neuman-Dirichlet 
type  perform  quite  poorly  for  certain  problems  of  the  form  (3);  cf.  Greenbaum,  Li 
and  Chao,  [30].  For  other  algorithms  a  satisfactory  theory  already  exists  with  bounds 
which  only  depend  on  the  local  variation  of  p(i);  see  Section  5.  For  stiU  others,  this 
issue  is  still  to  be  fully  addressed. 
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The  triangulation  of  ft  is  introduced  in  the  following  way.  The  region  is  first 
divided  into  non-overlapping  substructures  J),  ,  i  =  l,---,iV.  We  assume  that  the 
substructures  are  chosen  so  that  the  discontinuities  of  p(i)  occur  only  at  substructure 
boundaries.  To  simplify  the  description,  we  confine  most  of  our  study  to  triangular 
(simplicial)  substructures.  In  such  a  case,  the  original  region  must,  of  course,  be  a 
polygon  (polyhedron).  We  note  that  the  rate  of  convergence  of  the  iterative  methods, 
which  are  considered  in  this  paper,  depend  only  very  mildly  on  the  size  of  the  sub- 
structures. The  size  of  the  subregions  can  therefore  be  selected  primarily  to  minimize 
the  cost  of  a  single  iteration  on  a  particular  computer  system. 

All  the  substructures  Q,  are  further  divided  into  elements.  The  common  assump- 
tion in  finite  element  theory  that  all  elements  are  shape  regular  is  adopted  and  the 
same  assumption  is  made  concerning  the  substructures.  On  the  element  level  this 
means  that  there  is  a  uniform  bound  on  hfc/rn  ,  which  is  independent  of  the  number 
of  degrees  of  freedom.  Here  h[^  is  the  diameter  of  the  element  K  and  rj^-  the  diameter 
of  the  largest  inscribed  sphere  in  A'. 

Since  a(uh,Vh)  =  a(u,Vh),  V  u/,  G  V''',  the  finite  element  solution  is  the  projection 
of  the  exact  solution  onto  the  finite  element  space  with  respect  to  the  inner  product 
defined  by  the  bilinear  form.  The  problems  defined  on  the  subregions,  from  which  pre- 
conditioners  for  the  entire  problem  can  be  assembled,  can  often  similarly  be  viewed  in 
terms  of  orthogonal  projections  onto  subspaces  directly  associated  with  the  subregion 
in  question. 

We  also  need  to  use  matrix  representations  of  the  finite  element  problems  since 
many  algorithmic  details  can  best  be  described  using  matrix  notations.  When  doing 
this,  it  is  convenient  to  consider  domain  decomposition  in  light  of  structural  engi- 
neering computational  practices.    The  elements  of  the  stiffness  matrix  K  are  given 

by 

where  t^/  and  9^  are  standard  finite  element  basis  functions.  Since  an  integral  over  Q 
can  be  written  as  a  sum  of  integrals  over  the  substructures,  the  stiffness  matrix  can 
be  assembled  from  the  stiffness  matrices  7v'''',  which  have  the  elements 

4m  =  an,(V(,Vm)  • 

Here  the  relevant  /  and  m  correspond  to  degrees  of  freedom  a.ssociated  with  the  closure 
of  the  substructure  n,.  The  bilinear  form  aQ^{ufi,Vh)  represents  the  contribution  to 
the  integral  aQ{uh,Vh)  from  the  substructure  f2,. 

This  so-called  subassembly  process  can  be  summarized  in  the  formula 

(4)  i^A'2/  =  ^i(''^A'<')y('), 

where  x^'^  is  the  subvector  of  parameter  values  associated  with  H,. 

Remark.  In  practice,  the  subassembly  process  is  often  used  recursively,  creating 
larger  and  larger  so-called  super  elements.  This  process  can  be  interleaved  with  the 
elimination  of  all  the  variables  which,  at  the  level  in  question,  are  coupled  only  to 
other  variables  of  the  same  super  element;  cf.  e.g.  [4].  Except  for  a  brief  discussion  in 
Section  4,  we  consider  only  three  levels  in  this  paper:  the  elements  with  a  characteristic 
diameter  h,  the  substructures  with  a  diameter  on  the  order  of  H  and  the  entire  region 
fi,  which,  without  loss  of  generality,  is  assumed  to  have  unit  diameter. 
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If  we  divide  the  subvectors  i''^  cissociated  with  the  i-th  substructure  into  two,  Xj 
and  Xg  ,  corresponding  to  the  variables  which  are  interior  to  the  substructure  and 
those  which  are  shared  with  other  substructures,  then  the  matrix  A"'"'  can  be  written 


II        ^^IB 


^^IB        ^^BB 

Since  the  interior  variables  are  associated  with  only  one  of  the  substructures,  they 
can  be  eliminated  locally  and  in  parallel.  The  reduced  matrix  is  a  so-called  Schur 
complement  and  has  the  form 

•^        -  -'^BB       -'^/B    -'^//       ^^  IB- 
It  is  now  easy  to  show  that  if  the  corresponding  Schur  complement  of  the  global 
stiffness  matrix  K  is  denoted  by  5,  then, 

(5)  xlSyB  =  Y.^^fs(^^y^S- 

The  elimination  of  the  interior  variables  from  the  substructures  can  be  viewed  in 
terms  of  orthogonal  projections,  with  respect  to  the  bilinear  form,  of  the  solution  Uh. 
of  equation  (1)  onto  the  subspaces  nQ{Q.t){^V^,  i  =  l,---,N.  It  is  easy  to  show  that 
these  subspaces  are  orthogonal,  in  the  sense  of  the  bilinear  form  (4),  to  the  so-called 
piecewise  discrete  harmonic  functions  given  by 

A'gx^' +  A'i']xg  =  0,    \fi. 

If  the  local  problems  are  solved  exactly,  what  remains  is  to  find  a  sufficiently  ac- 
curate approximation  of  the  part  of  the  solution  which  is  piecewise  discrete  harmonic. 
This  can  be  done  by  approximately  solving  the  reduced  Linear  system  with  the  matrix 
5.  A  specific  iterative  substructuring  method  is  obtained  by  selecting  a  precondi- 
tioner  for  the  matrix  5.  Once  an  approximation  of  the  solution  has  been  found  on  the 
boundaries  of  the  substructures,  the  solution  can  be  found  everywhere  by  separately 
solving  local  Dirichlet  problems  on  each  substructure. 

An  important  family  of  domain  decomposition  algorithms  can  be  derived  by  re- 
placing each  contribution  5*'^  to  the  Schur  complement  by  a  different  matrix  5*''.  A 
preconditioner  5  of  5  is  then  created  by  assembling  the  local  contributions  in  the  same 
way  as  in  equation  (5).  When  choosing  S,  it  is  important  to  create  a  linear  system, 
which  is  cheaper  to  solve  and,  at  the  same  time,  to  ensure  that  k{S~^S)  is  small.  A 
natural  idea  is  to  eliminate  the  coupling  between  the  groups  of  variables  associated 
with  different  edges  (faces)  of  the  individual  substructures.  Further  details  are  given 
in  Dryja  and  Widlund  [28],  where  this  process  is  interpreted  as  a  splitting  in  the  sense 
of  Varga  [51].  In  the  analysis,  a  simple  and  powerful  idea  is  based  on  the  fact  that  if 

(7)  c.x(^)^5(")x(^)    <  x^fs^'^x^^    <  Qx^fs^'^x^S  ,  Vx^', 
then, 

(8)  cxgSxB    <x'^SxB   <Cx%SxB,yxB, 


U"+l/2 

=     / 
=     u" 

in  O'l, 

-Au"+i     = 

/ 

lt"+l/2 

in  n'2, 
ondfl'2 

where 

(9)  c    =    min,  c,     and     C    =    max,  C,  . 

When  this  technique  can  be  used,  it  is  as  easy  to  derive  satisfactory  bounds  for 
equations  (3)  as  for  (2);  cf.  Section  5. 

3.  Multiplicative  and  Additive  Schwarz  Methods.  In  this  section,  we  first 
review  Schwarz's  classical  alternating  method  [47],  its  variational  formulation,  cf.  P.  L. 
Lions  [34],  and  an  additive  variant  of  the  algorithm,  cf.  Dryja  [21],  Dryja  and  Widlund 
[26]  and  Matsokin  and  Nepomnyaschikh  [41].  We  then  discuss  a  general  approach  to 
the  analysis  of  the  rate  of  convergence  of  the  additive  algorithms. 

We  begin  by  briefly  discussing  the  classical  formulation  of  Schwarz's  method  in  the 
case  of  the  continuous  Poisson  equation.  There  are  two  fractional  steps  corresponding 
to  two  overlapping  subregions,  Q[  and  ^'2,  the  union  of  which  is  the  region  fi  .  Let  an 
initial  guess  u°  G  V  be  given.  The  iterate  u"''"^  is  determined  from  u"  by  sequentially 
updating  the  approximate  solution  in  the  two  subregions: 


and 


We  could  just  as  well  have  written  down  the  finite  element  version  of  the  algorithm. 
From  now  on,  we  only  consider  that  case.  It  is  easy  and  convenient  to  describe  this 
method  in  terms  of  two  projections  P,,  i  =  1,2,  onto  V/^  =  nQ{il'i)C\V^  ;  cf.  Lions 
[34].  The  projections  are  defined  by 

(10)  a{P,Vh,M  =  a(vH,4>h),  y4>h  e  Vt  . 

It  is  also  easy  to  show  that  the  error  propagation  operator  of  this  multiplicative 
Schwarz  method  is 

[I  -  P2){I  -  P,). 

This  algorithm  can  therefore  be  viewed  as  a  simple  iterative  method  for  solving 

(11)  (Px  +  P2  -  P2P^)Uh  =  gK, 

with  an  appropriate  right-hand  side  gh  ■ 

This  algorithm  can  be  generalized  immediately  to  any  number  of  subspaces.  Let 

V^  =  V^  +  V2  +  ...+  Vn  . 

The  subspaces  V^  can  be  chosen  quite  arbitrarily.   The  projections  are  defined  as  in 
equation  (10). 

The  operator  of  (11)  is  a  polynomial  of  degree  two  and,  thus,  the  method  is  not 
ideal  for  parallel  computing,  since  two  sequential  steps  are  involved.  If  more  than 
two  subspaces  are  used,  this  effect  is  further  pronounced  even  if  the  degree  of  the 
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polynomial  representing  the  multiplicative  algorithm  often  is  lower  than  maximal. 
This  is  so  because  a  product  of  two  projections  associated  with  subspaces  that  do 
not  intersect  nontrivially  vanishes;  cf.  further  Widlund  [52],  [55].  This  can  best 
be  described  in  terms  of  the  coloring  of  an  undirected  graph  in  which  the  nodes 
represent  the  subspaces  and  the  edges  nontrivial  intersections  of  pairs  of  subspaces. 
The  problems  associated  with  all  subspaces  of  the  same  color  can  be  solved  in  parallel 
since  those  subspaces  are  mutually  orthogonal.  For  the  purpose  of  theory,  we  can  then 
regard  all  subspaces  of  the  same  color  as  one  subspace;  this  can  greatly  improve  the 
upper  bound  for  the  eigenvalues  of  the  operator  P,  which  we  are  about  to  introduce. 
The  basic  idea  behind  the  additive  form  of  the  algorithm  is  to  work  with  the 
simplest  possible  polynomial  in  the  projections.  The  equation 

(12)  Puh  =  {Pi  +  P2-V....  +  PN)uh  =  9'h, 

is  solved  by  an  iterative  method.  (In  fact,  we  could  just  as  well  work  with  Yl^=\  c^t-fn 
where  a,  >  0  are  suitably  chosen  constants.)  Since  we  can  show  that  the  operator  P 
is  symmetric  and  positive  definite,  with  respect  to  the  bilinear  form,  the  method  of 
choice  is  the  conjugate  gradient  method.  Equation  (12)  must  have  the  same  solution  as 
equation  (1),  i.e.  the  correct  right-hand  side  must  be  found.  Since  a(uti,4>h)  =  fi4>h)  i 
by  equation  (1),  the  right-hand  side  g'f^  can  be  constructed  by  solving  equation  (10)  for 
all  values  of  i  and  adding  the  results.  It  is  similarly  possible  to  apply  the  operator  P  of 
equation  (12)  to  any  element  of  V^  by  applying  each  projection  P,  to  the  element  and 
adding  the  results.  Most  of  the  work,  in  particular  that  which  involves  the  individual 
projections,  can  be  carried  out  in  parallel.  We  note  that  an  additive  Schwarz  algorithm 
is  fuJly  specified  by  its  subspaces.  In  this  survey  paper,  we  will  sometimes  just  specify 
the  subspaces  without  much  discussion  of  how  the  corresponding  algorithm  can  be 
implemented. 

It  is  well  known  that  the  number  of  steps  required  to  decrease  an  appropriate 
norm  of  the  error  of  a  conjugate  gradient  iteration  by  a  fixed  factor  is  proportional  to 
y/n  ,  where  k  is  the  condition  number  of  the  relevant  operator;  see  e.g.  Golub  and  Van 
Loan  [29]  .  We  therefore  need  to  establish  that  the  operator  P  of  equation  (12)  is  not 
only  invertible  but  also  that  satisfactory  upper  and  lower  bounds  on  its  eigenvalues 
can  be  obtained. 

An  upper  bound  for  the  eigenvalues  of  P  is  given  by  N  since  P  is  the  sum  of 
projections.  By  combining  all  subspaces  of  the  same  color  into  one,  as  indicated 
above,  this  bound  can  be  improved  by  replacing  A'^  by  the  number  of  colors  of  any 
coloring  of  that  graph.  For  several  algorithms  this  results  in  an  upper  bound,  which 
is  independent  of  the  number  of  subspaces.  In  other  instances,  a  useful  technique 
is  provided  by  strengthened  Cauchy  inequalities;  cf.  Mandel  and  McCormick  [38], 
Widlund  [56]  and  Yserentant  [58]. 

A  lower  bound  can  often  be  obtained  conveniently  by  using  a  lemma,  inspired  by 
Lions  [34];  the  simple  proof  is  also  given  in  Widlund  [56]. 

Lemma  1.  Let  Uh  =  Xltli  "'i.i'  where  Uh^i  e  V,,  be  a  representation  of  an  element 
ofV^  =  Vi  +  ...  +  Vn.  If  the  representation  can  be  chosen  so  that  ^Z^i  a.{uh  ,,  Uh  ,)  < 
Cla{uH,Uh),'iuh  eV^,  then  Xrmr^i P)>C^'^. 

Remark.  If  we  expand  individual  subspaces,  there  is  a  larger  choice  in  selecting 
Uh,i  6  Vf  The  best  bound  for  Co  can  then  only  improve.  If  we  can  expand  the 
subspaces  without  worsening  the  upper  bound,  which  is  often  possible,  our  estimate 
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of  k(P)  improves.  On  the  other  hand,  an  enlarged  subspace  also  means  that  the 
subproblem  has  more  variables  and  that  it  is  worse  conditioned.  For  the  special  case 
of  the  classical  Schwarz  method  on  two  regions,  this  tradeoff  is  well  understood;  for  a 
discussion  of  precise  estimates  of  the  rate  of  convergence;  see  Bjorstad  and  Widlund 
[10]. 

Remark.  In  his  recent  dissertation.  Tarek  Mathew  [40]  has  given  a  bound  on  the 
rate  of  convergence  of  the  multiplicative  algorithm  in  terms  of  Co  and  the  number 
of  colors  of  the  graph.  While  his  bound  can  probably  be  improved  considerably,  it 
nevertheless  shows  that  if  these  parameters  are  bounded  independently  of  the  mesh 
parameters,  then  the  spectral  radius  of  the  error  propagation  operator  is  uniformly 
bounded  by  a  constant  that  is  less  than  1. 

One  of  the  attractive  features  of  the  framework  introduced  in  this  section  is  the 
ease  by  which  variants  of  the  basic  algorithm,  obtained  by  replacing  the  local  linear 
systems  by  preconditioners,  can  be  analyzed.  Let  us,  for  example,  consider  the  additive 
Schwarz  method  for  a  finite  element  approximation  of  the  problem  discussed  in  the 
beginning  of  this  section.  We  can  write  the  projection  Pi  in  matrix  terms.  After  a 
suitable  permutation  of  the  variables,  it  corresponds  to 

(13)  y  =  P,x=(^^'f)     °)a-x. 

Here  A'(i)  represents  the  stiffness  matrix  of  the  Dirichlet  problem  on  Cl[.  It  is  easy  to 
see  that  the  matrix  of  (13)  is  symmetric  in  the  /iT-inner  product  that  corresponds  to 
the  bilinear  form.  If  A'j~|  is  replaced  by  A'j~>,  etc.,  then  it  is  easy  to  show  by  using  a 
Rayleigh  quotients  argument  that  the  eigenvalues  of  the  resulting  operator  P  satisfy 

(14)  A„^n(P)mm(A^„(A'(.)A',-J))    <    A„^„(P), 
and 

(15)  Amax(^)    <    >^Tn^x{P)raa.x{X^^{Ki^,^K-^))  . 

An  estimate  of  k(P)  follows  immediately, 

max,  Amax(A'(,)A'r^^) 

(16)  K(P)<   !_LJil_K(P). 

mm,  Amin(A'(i)A'(,)  ) 

4.  Additive  Schwarz  Methods  on  Overlapping  Subregions.  We  now  de- 
scribe the  additive  Schwarz  method  introduced  in  Dryja  and  Widlund  [26],  [28];  cf. 
also  Dryja  [21].  In  those  papers  only  Poisson's  equation  was  considered.  In  the  origi- 
nal technical  report,  good  results  were  given  only  for  problems  in  two  dimensions,  but 
we  now  need  no  such  restrictions.  The  most  satisfactory  proof,  so  far,  is  given  in  [28]. 
Here,  we  also  discuss  problems  of  the  form  (3). 

We  start  with  the  same  triangular  (simplicial)  non-overlapping  substructures  ft,- 
that  have  been  considered  before.  We  extend  each  substructure  to  a  larger  region 
ft'We  assume  that  the  overlap  is  generous,  assuming  that  the  distance  between  the 
boundaries  5f),  and  d^'^  is  bounded  from  below  by  a  fixed  fraction  of  Hi,  the  diameter 
of  ft,.  We  also  assume  that  dil[  does  not  cut  through  any  element.  We  make  the  same 
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construction  for  the  substructures  that  are  next  to  the  boundary  except  that  we  cut 
off  the  part  of  0,'-  that  is  outside  of  Q  . 

Remark.  The  analysis  of  Schwarz  methods  is  more  complicated  when  the  overlap 
is  less  generous;  cf.  a  discussion  in  Lions  [34].  Such  a  situation  occurs  if  the  region 
is  L-shaped  and  partitioned  into  two  overlapping  rectangles.  A  similar  situation  is 
discussed  in  Dryja  and  Widlund  [28]  in  an  analysis  of  iterative  substructuring  methods; 
cf.  also  Section  5. 

Our  finite  element  space  is  represented  as  the  sum  of  N  +  1  subspaces 

v^  =  Vo^  +  vl"  + ...  +  y/} . 

The  first  subspace  Vq*  is  equal  to  V^,  the  space  of  continuous,  piecewise  linear  func- 
tions on  the  coarse  mesh  defined  by  the  substructures  fi,  .  The  other  subspaces  are 
related  to  the  subdomains  in  the  same  manner  as  in  a  traditional  Schwarz  algorithm, 
i.e.  V>^ -^  n^iil[)f]V^  . 

Remark.  There  are  few  differences  between  the  problem  related  to  the  first  sub- 
space  and  the  others.  The  right  hand  side  of  the  particular  linear  system  is  generated 
as  weighted  averages  with  weights  determined  by  the  coarse  mesh  basis  functions.  We 
also  need  to  interpolate  the  solutions  in  Vq*  when  combining  them  with  the  contri- 
butions from  the  other  subspaces.  The  coarse,  global  approximation  of  the  elliptic 
equation  is  otherwise  quite  similar  to  the  local  problems. 

The  global  coarse  problem  provides  a  mechanism  for  the  global  transportation 
of  information.  As  shown  in  Widlund  [54],  the  rate  of  convergence  of  any  domain 
decomposition  method,  for  which  the  interaction  is  only  through  next  neighboring 
subdomains,  has  a  condition  number  which  grows  at  least  as  fast  as  1/ H^. 

The  following  result  is  established  in  Dryja  and  Widlund  [28]. 

Theorem  1.  The  operator  P  of  the  additive  algorithm  defined  by  the  spaces  V^ 
and  V,  ,  applied  to  equation  (2),  satisfies  the  estimate  k(P)  <  const. 

In  the  proof,  a  quasi-interpolant,  ///  ,  which  is  a  bounded  operator  in  ^o(^)  ^^'^ 
satifies 

\\uh  -  iHUhWL^iQ)  <  const. H\uh\f{^n), 

is  used.  By  examining  the  proof,  it  is  easy  to  see  that  if  uniform  bounds  of  the  same 
form  hold  for  the  weighted  spaces  L2,pi^)  and  HopiQ),  defined  by 

/    p{x)\u\'^dx      and        /  p(x)\Vu\'^dx, 

respectively,  then  the  optimality  also  holds  for  the  variable  coefficient  case  (3).  This  is 
a  difficult  problem,  still  not  fully  understood;  cf.  Xu  [57]  for  partial  results.  We  note 
that  for  two-dimensional  problems  the  proof  given  in  Dryja  and  Widlund  [26]  works 
just  as  well  for  equation  (3).  The  bound  of  the  condition  number  given  in  that  paper 
is  proportional  to  (1  +  logi H/h)). 

We  have  also  considered  the  use  of  more  than  two  levels  of  partitioning  of  the  re- 
gion and  established  that  for  equation  (2)  k(P)  does  not  grow  faster  than  quadratically 
in  the  number  of  levels.  We  do  not  know  if  this  result  can  be  improved.  For  convex 
regions,  we  have  shown  that  the  growth  is  linear  in  the  number  of  levels.  This  result 
cannot  be  improved.  Our  work  has  been  carried  out  jointly  with  Xuejun  Zhang,  who 

9 


has  also  made  a  large  number  of  numerical  experiments.  Rapid  convergence  can  be 
obtained  when  using  such  methods  solving  only  very  small  linear  systems  of  equations. 
As  we  have  already  pointed  out  in  the  introduction,  Xiao-Chuan  Cai  [20],  [19]  and 
Tarek  Matew  [40]  have  extended  our  analysis  to  certain  nonsymmetric  and  indefinite 
problems. 

5.  Iterative  Substructuring  Methods  as  Schwarz  Methods.  In  this  sec- 
tion, we  show  how  iterative  substructuring  methods  can  be  derived  by  using  the 
framework  of  subspaces  and  projections  developed  in  Section  3.  We  first  consider 
problems  in  the  plane  and  give  an  estimate  of  the  condition  number  of  the  operator 
P  that  corresponds  to  a  specific  choice  of  subspaces.  For  this  basic  algorithm,  the 
proof  of  the  estimate  does  not  rely  on  the  extension  theorems,  which  have  played  an 
important  role  in  previous  estimates  for  iterative  substructuring  methods;  cf.  Wid- 
lund  [53].  We  then  indicate  how  results  for  the  special  case  of  two  substructures,  and 
an  argument  about  preconditioners  for  the  subspace  problems,  can  be  used  to  derive 
different  algorithms,  including  one  due  to  Bramble,  Pasciak  and  Schatz  [15].  We  also 
consider  a  second  choice  of  subspaces  and  show  that  the  resulting  algorithm  extends 
to  a  fast  method  for  the  three-dimensional  czise.  The  resulting  algorithm  is  similar  but 
not  identical  to  one  of  the  two  algorithms  considered  in  Bramble,  Pasciak  and  Schatz 
[17];  cf.  also  Dryja  [24]  for  early  work  on  the  three-dimensional  ca^e.  The  discussion 
in  this  section  extends  that  in  Dryja  and  Widlund  [28]  and  relies  in  part  on  ideas  from 
Mandel  [36]. 

We  assume  that  the  region  is  divided  into  substructures  and  elements  as  in  Section 
2.  We  first  consider  the  two-dimensional  case  and  introduce  an  additive  Schwarz 
method.  We  use  the  coarse  space  V^  introduced  in  Section  4  and  let  the  subregions 
f^ij  =  ^t  U  ^t]  U  ^j  pls-y  ^he  same  role  as  Q.\  in  Section  4.  Here  fi,  and  Vij  are  adjacent 
substructures  with  a  common  edge  ?,_,.  The  local  subspaces  V^^  are  thus  defined  as 

Compared  with  the  subspaces  used  in  the  previous  section,  we  use  less  overlap 
in  the  sense  that  only  the  elements  of  V^  can  differ  from  zero  at  the  vertices  of  the 
substructures.  This  is  reflected  in  a  poorer  bound  on  the  condition  number  when  we 
work  with  Lemma  1. 

Theorem  2.  In  the  two-dimensional  case,  the  operator  P  of  the  additive  algorithm 
defined  by  the  spaces  V^  and  V^,j  satisfies  the  estimate  k{P)  <  const.{l  +  log{n/h))^ 
for  equations  (2)  and  (3).  The  constant  is  independent  of  h,  H  and  the  discontinuities 
of  the  coefficient  p{x). 

For  this  algorithm,  Xm&xiP)  <  4  ;  cf.  Dryja  and  Widlund  [28].  It  is  ecisy  to  see 
that  no  more  than  four  colors  are  needed  to  color  the  graph;  each  interior  substructure 
is  covered  exactly  three  times  by  the  subspaces  associated  with  its  three  edges  and, 
in  addition,  we  have  the  global  space  V^  which  intersects  all  other  subspaces.  In  our 
proof  of  the  lower  bound  of  the  spectrum  of  P,  given  fully  in  [28],  we  use  Lemma  1 
and  the  following  lemma,  which  also  plays  an  important  role  in  the  more  traditional 
theory  of  iterative  substructuring  algorithms. 

Lemma  2.   Let  a  be  any  convex  combination  of  values  of  Uh{x),x  G  fi,.  Then 

\Wh  -  a|lioo(n,)  <  const.  [\  -t-  log(i7//i))|uh|^,^Q_j  . 
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Variations  of  this  result  date  back  at  least  to  1966  and  proofs  are  given  in  a  number 
of  papers;  see  e.g.  Bramble  [12],  Bramble,  Pcisciak  and  Schatz  [15]  or  Yserentant  [58]. 

Since  all  the  elements  of  V,^  vanish  at  the  vertices  of  the  substructures,  we  must 
choose  the  interpolant  ///u/i  as  the  element  in  V^  when,  as  required  by  Lemma  1,  we 
define  the  representation  of  u/, .  It  is  easy  to  show  that  I^z/'^^I^wQ.n  can  be  estimated 
from  above  by 

(17)  ^  (u;,(IO-w/t(a:/))^ 

where  the  V^  is  the  set  of  vertices  of  the  substructure  Cl, .  It  then  follows  from  Lemma 
2  that 

(18)  \^fi''>^\]mQ,)  ^  con5i.(l  +  login/ h))\uk\]^,^^^^  . 

For  a  bound  on  the  other  terms  in  the  decomposition,  required  when  using  Lemma  1, 
see  Dryja  and  Widlund  [28]. 

Using  the  method  with  local  preconditioners  developed  in  the  end  of  Section  3,  a 
number  of  algorithms  can  be  derived  from  the  basic  method  that  we  have  just  intro- 
duced. All  that  is  needed  is  to  replace  the  Dirichlet  problem  on  ri,j  by  a  preconditioner 
with  a  symmetric,  positive  definite  coefficient  matrix.  A  proof  of  the  main  result  of 
Bramble,  Pasciak  and  Schatz  [15]  can  be  obtained  in  this  way  by  using  (14),  (15), 
Theorem  2  and  a  bound  for  the  condition  number  of  the  two  subregion  problem  de- 
fined on  ilij.  Such  local  bounds  are  given  e.g.  in  Bjorstad  and  Widlund  [9],  Bramble, 
Pasciak  and  Schatz  [18]  and  Dryja  [22],  [23]. 

In  our  previous  paper,  only  equation  (2)  was  considered.  By  using  the  formulas 
(7),  (8),  and  (9),  the  proof  can  be  extended  immediately  to  equation  (3);  the  whole  ma- 
chinery can  also  be  used  locally.  In  this  variant  of  the  proof,  only  four  subspaces  play  a 
role,  namely,  the  three-dimensional  space  of  linear  functions  on  fi,,  and  the  restriction 
to  Qi  of  the  spaces  V-^  which  correspond  to  the  three  edges  of  that  substructure. 

In  preparation  for  the  three-dimensional  case,  we  develop  an  alternative  basic 
algorithm  for  the  two-dimensional  problems.  This  algorithm  is  defined  in  terms  of 
subspaces  related  to  individual  substructures.  For  any  substructure  fi,,  we  keep  the 
three  local  subspaces  obtained  as  above  from  the  spaces  F,^.  For  an  interior  substruc- 
ture, we  replace  the  three-dimensional  space  of  linear  functions,  used  previously,  by 
the  three-dimensional  space  spanned  by 

( 19)  4'^  =  ^k-  tf '  E  V'  +  tf ' '  f'  e  v.. 

'ev, 


Here  ipk  is  the  standard  nodal  basis  function  associated  with  an  element  of  V,-.  We  note 
that,  just  as  the  regular  nodal  basis  functio: 
constants  /3^    are  chosen  below,  so  that 


(•^^-^  =  Ski,    k,l  6  v..  The  positive 


kev, 
Then, 


kev, 


E  4''  =  1. 

JtGV, 
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This  subspace,  thus,  also  contains  the  constants,  i.e.  the  null  space  of  the  elliptic 
operator  restricted  to  Q,  with  a  Neumann  boundary  condition.  (When  these  methods 
are  extended  to  elasticity  problems,  we  similarly  have  to  include  the  whole  null  space  of 
that  operator  in  the  space  that  accounts  for  the  global  transportation  of  information.) 
For  a  substructure  that  intersects  the  boundary,  we  simply  use  the  space  spanned 
by  the  standard  basis  functions  ^k  corresponding  to  the  vertices  of  the  substructures 
where  a  Dirichlet  condition  is  not  imposed.  For  the  correct  definition  of  boundary 
substructures  for  the  case  of  three  dimensions;  see  below. 

An  elementary  computation  shows  that,  for  an  interior  substructure,  the  quadratic 
form  associated  with  this  subspace  is  given  by 

(20)  E4'2("'.(^fc)-4*'')'- 

Here,  ir^  is  the  weighted  average  of  Uh{xk)-,  k  £  V,-,  with  the  weights  ^[''.  It  is  easy  to 
see  that  this  quadratic  form  can  be  bounded  from  above  and  below  by  the  form  (17) 
given  in  the  discussion  of  the  previous  algorithm. 

Inspired  by  Mandel  [36],  we  now  introduce  the  quadratic  form 


(21)  J{uh,y)  =  E  L  ^^^kli^hixk)  -  t/<*')'  -  2l^  E  ^hixk)b, 


ca(u,(x,)-t/('))2-2x:^^'.^^.^A(') 

kev,  1    kev, 


For  substructures  which  have  at  leeist  one  vertex  on  dCl,  we  set  y^'^  =  0.  In  three 
dimensions,  we  adopt  the  same  rule  for  all  substructures  which  have  at  least  an  edge 
on  dCl,  but  treat  those  which  only  touch  the  boundary  at  individual  points  as  if  they 
were  interior. 

The  vector  6(''  corresponds  to  the  part  of  the  load  vector  which  is  associated  with 
n,.  Its  components  are  computed  by  evaluating  an  appropriate  functional  at  tp)^  . 

It  is  easy  to  show  that  if 

^^k     "  T7\   ' 

then, 

(22)  4')  =  argmm  ^  42(t.;.(i^)  -  y^'^)\ 

Therefore,  the  solution  of  the  linear  system  corresponding  to  these  subspaces  can  be 
obtained  by  minimizing  J  with  respect  to  Uh  and  y.  By  setting  the  gradient  of  J 
equal  to  zero,  we  obtain  a  linear  system  of  equations,  which  naturally  can  be  written 
in  a  two-by-two  block  form.  The  diagonal  blocks  are  diagonal  with  elements  IZi^il 
^•nd  ^A-gv  kkk'  respectively.  The  off-diagonal  blocks  are  sparse  and  are  also  given  by 

'^kk- 

Since  the  diagonal  blocks  are  diagonal,  it  is  easy  to  eliminate  all  the  variables 
f^hixk)-  This  elimination  step  results  in  a  sparse  system,  similar  to  that  of  a  finite 
difference  problem  on  a  coarse  mesh  that  is  dual  to  the  one  used  to  define  the  space 
V" .  It  is  often  feasible  and  economical  to  solve  this  system  by  a  direct  method.  Once 
the  values  ofu]^  are  known,  it  is  easy  to  determine  the  values  of  u/,(ijt)>^  G  V,,  and  the 
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corresponding  element  of  the  subspace.  The  contributions  of  the  other  subspaces  can 
now  be  computed,  while  observing  that  we  must  make  the  sum  of  all  the  contributions 
continuous  across  the  interfaces  between  the  substructures. 

Remark.  The  method  just  outlined  draws  heavily  on  the  work  of  Bramble,  Pasciak 
and  Schatz  [16].  However,  the  current  algorithm  uses  a  different  average  and  basis. 
Our  theory  does  not  require  that  the  linear  system,  which  determines  the  averages,  to 
have  a  M  — matrix.  Our  method  can  therefore  be  extended  to  the  elliptic  systems  of 
elasticity;  cf.  also  Mandel  [36]. 

It  can  be  shown,  straightforwardly,  that  there  is  no  longer  a  need  to  insist  on 
triangular  substructures.  In  the  general  case,  we  can  simply  redefine  V,  as  the  set  of 
nodal  points  which  belong  to  the  boundaries  of  at  least  three  substructures. 

Lemma  2  holds  only  in  two  dimensions;  it  resembles  a  Sobolev  inequality  which 
is  far  from  valid  in  three  dimensions.  However,  for  finite  element  functions,  we  can 
find  a  useful  bound  of  the  L2  -  norm  over  an  interval  in  terms  of  the  strain  energy. 
Let  W,  denote  the  wire  basket  of  the  substructure  H,.  This  is  the  union  of  the  edges 
of  the  tetrahedral  substructure;  more  generally,  it  is  defined  as  the  set  of  nodal  points 
of  fi,,  which  belong  to  the  boundaries  of  at  least  three  substructures.  The  following 
lemma  is  essentially  a  corollary  of  Lemma  2. 

Lemma  3.   Let  a  be  any  convex  combination  of  values  of  Uh(x),x  £  W,.  Then 

h\\uh  -  ajlfjjyv^)  <  const.  (1  +  log{H/h))\uh\]ji^Q^^  ■ 


An  algorithm  for  the  three-dimensional  case  can  now  be  defined  in  terms  of  the 
subspaces  V/j,  which  are  defined  just  as  in  the  two-dimensional  case,  and  a  special 
subspace.  This  subspace  is  given  in  terms  of  a  basis  defined  as  in  (19),  where  we 
replace  the  set  V,  by  W,.  The  study  of  the  linear  algebra  problem  associated  with 
this  basis  proceeds  very  much  as  in  the  two-dimensional  case.  We  only  note  that 
when  we  compute  the  expression  corresponding  to  (20),  we  also  get  off-diagonal  terms 
corresponding  to  nodes  on  VV^,  which  are  next  neighbors.  It  is  however  easy  to  show 
that  these  contributions  can  be  disregarded.  The  rest  of  the  analysis  and  the  devel- 
opment now  follows  just  as  in  the  two-dimensional  case.  We  note  that  the  number  of 
elements  in  W,  grows  linearly  with  H/h.  The  nontrivial  part  of  the  sparse  linear  sys- 
tem of  the  special  subspace  has  a  dimension,  which  is  equal  to  the  number  of  interior 
substructures. 

We  have  obtained  the  following  result. 

Theorem  3.  The  operator  P  of  the  additive  algorithm  defined  by  the  spaces  V/^ 
and  that  defined  by  the  basis  functions  given  by  (19),  or  the  corresponding  formula  for 
the  three-dimensional  case,  satisfies  the  estimate  k{P)  <  const.{l  +  log(H/h))^  for 
equations  (2)  and  (S).  The  constant  is  independent  of  h,  H  and  the  discontinuities  of 
the  coefficient  p{x). 

6.  Schwarz  Methods  for  Problems  on  the  Interfaces.  We  have  already  seen 
that  the  original  linear  systems  of  equations  can  be  reduced  to  smaller  systems,  where 
the  remaining  variables  are  associated  with  the  boundaries  89.^  of  the  substructures. 
We  have  also  shown  that  this  reduced  set  of  variables  are  associated  with  the  piece- 
wise  discrete  harmonic  functions.  This  is  quite  similar  to  potential  theory  for  elliptic 
problems  where  there  is  also  a  reduction  in  dimension. 
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Here,  we  only  consider  the  first  model  case  (2),  since  at  this  time,  we  do  not 
know  if  bounds  independent  of  the  variations  of  the  coefficients  of  equation  (3)  can 
be  obtained  for  the  methods  that  we  are  going  to  introduce.  We  denote  by  V^^^^{T) 
the  space  of  piecewise  discrete  harmonic  functions  defined  by  the  values  of  the  finite 
element  functions  on  the  set 

(23)  T  ^\Jdn,\dn. 

The  space  of  traces  H^/'^{T)  are  the  restrictions  of  n^{?l)  to  F.  As  is  well  known,  the 
norm  of  this  space  can  be  defined  by 

(24)  l"l//i/2(r)  =  Z^l"l//i/2{9n.)' 
where 

(25)  W\mnidu,)=  I      I     {Hx{s))-u{x{t))\^l\x{s)-x{t)\')dsdt. 

Here,  J  =  2,3,  for  problems  with  fi  C  R"^. 

We  consider  several  different  additive  Schwarz  methods.  The  first  is  defined  by 
the  partitioning 

(26)  yLmi^)  =  Y.^Lm,i. 

i 

where  V^''^^^  ,  is  the  subspace  of  ^';i''arm(^)  "^W-h.  zero  values  on  F  outside  5f2,.  By  using 
the  same  technique  as  in  the  previous  three  sections,  it  is  easy  to  see  that  the  operator 
P,  which  correspond  to  this  set  of  subspaces,  is  bounded  uniformly  from  above.  We  can 
use  the  techniques  of  Widlund  [54]  and  Dryja  [24]  for  the  two-  and  three-dimensional 
cases,  respectively  and  Lemma  1,  to  obtain  the  bound  Cq  <  const.{l  +  log{n /h))"^.  It 
is  not  possible  to  obtain  a  uniform  bound,  since  then  the  spaces  ^qo  ^^^  H^^^  would 
have  to  be  the  same;  cf.  Lions  and  Magenes  [33]. 

The  computations  associated  with  the  space  V/^^^^ ,  involve  piecewise  discrete 
harmonic  functions,  which  differ  from  zero  in  all  the  substructures  which  are  next 
neighbors  of  Q,.  The  coefficient  matrix  of  the  linear  system,  which  is  solved  when 
computing  the  projection  of  V,^^^^  onto  V^^^^-  is  the  principal  minor  of  S  cissociated 
with  the  nodes  on  dQ,.  It  also  involves  contributions  from  the  matrices  S^-'^  of  the 
neighbors  of  Q,  .  It  is  natural  to  replace  this  matrix  by  S^'\  since  the  use  of  this 
preconditioner  only  involves  a  problem  on  one  substructure. 

What  we  have  just  given  is  an  alternative  derivation  of  an  algorithm  introduced 
by  Bourgat,  Glowinski,  Le  Tallec  and  Vidrascu  [11].  We  note  that  there  is  a  minor 
technical  issue  related  to  the  fact  that  S^'^  is  singular  for  interior  substructures,  but 
that  it  is  easy  to  solve  that  problem.  Using  the  formalism  developed  in  Section  2,  we 
can  write  down  the  following  formula  for  the  preconditioner  which  corresponds  to  this 
algorithm, 

(27)  x^5-it/B  =  E^BV^"^i'^- 

Using  Lemma  3.2  in  Widlund  [54],  the  corresponding  estimates  developed  for  the 
three-dimensional  case  in  Dryja  [24]  and  the  estimate  (16),  we  can  conclude  the  proof 
of  the  following  result. 
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Theorem  4.  The  operator  P  of  the  additive  algorithm  defined  by  the  spaces 
Vharm,i,  applied  to  equation  (2),  and  the  local  preconditioner  defined  by  S^'^  satisfies 
the  estimate  k{P)  <  const. {I +  \og( H / h))'* .  The  constant  in  the  estimate  is  independent 
of  h  but  not  of  H. 

By  using  an  elementary  argument  given  in  Widlund  [54],  it  is  easy  to  see  that  k{P) 
must  grow  at  least  as  fast  as  1/ H^,  since  in  each  iteration,  information  is  exchanged 
only  with  neighboring  substructures.  The  performance  of  the  algorithm,  for  the  case 
of  many  substructures,  can  be  improved  by  introducing  the  same  spaces  of  modest 
dimension,  which  provide  some  global  transportation  of  information  in  the  methods 
considered  in  Section  5.  After  this  modification,  the  constant  in  the  theorem  can  be 
shown  to  be  independent  of  ^  as  well. 

We  note  that  in  each  iteration  of  this  algorithm,  both  a  Neumann  and  a  Dirichlet 
problem  have  to  be  solved  on  each  substructure.  This  results  in  twice  as  much  work 
per  iteration  as  for  many  other  domain  decomposition  methods. 

K  there  is  a  red-black  ordering  of  the  substructures,  there  is  another  way  of  an- 
alyzing this  algorithm.  We  first  consider  the  two  subdomain  case.  We  note  that  the 
inverse  of  the  preconditioner  corresponding  to  the  Neumann-Dirichlet  algorithm  is 
5'^'  or  S^^^  depending  on  the  roles  assigned  to  the  substructures;  cf.  Bj0rstad 
and  Widlund  [9].  We  can  therefore  view  the  preconditioner  introduced  in  this  section 
as  an  average  of  Neumann-Dirichlet  operators.  In  this  simple  case,  the  condition  num- 
ber remains  uniformly  bounded  since  there  are  no  so-called  cross  points.  Similarly,  we 
can  use  a  result  by  Dryja  [24]  to  prove  that  the  algorithm  for  three  dimensions,  which 
incorporates  the  global  subspace  introduced  in  the  previous  section,  has  a  condition 
number  which  is  bounded  by  const.{l  +  log(II/h))^.  This  again  follows  from  the  obser- 
vation that  this  new  preconditioner  can  be  viewed  as  an  average  of  two  instances  of  the 
operator  which  arises  in  the  study  of  the  version  of  the  Neumann-Dirichlet  algorithm 
considered  in  Dryja  [24]. 

We  conclude  this  section  by  a  brief  discussion  of  some  recent  work  by  Barry  Smith, 
a  Courant  Institute  graduate  student.  His  algorithm,  described  in  Smith  [48],  is  an 
additive  Schwarz  method,  which  uses  the  subspace  V^  both  in  two  and  three  dimen- 
sions. We  give  some  details  only  for  the  three-dimensional  case.  In  addition  to  V^, 
there  are  subspaces  associated  with  each  face,  edge  and  vertex  of  the  substructures. 
The  subspace  associated  with  a  face  is  .^qo  (rtj)  f]  ^harm  ^"^  therefore  closely  related 
to  the  space  F,^  introduced  in  Section  5.  Similarly,  for  an  edge  (vertex).  Smith  uses 
a  space  associated  with  a  neighborhood  in  F  of  the  edge  (vertex)  which  extends  a  dis- 
tance on  the  order  of  the  diameter  of  fi,.  The  success  of  this  approach  can  be  explained 
informally  by  noting  that  the  overlap  is  relatively  generous.  The  numerical  results, 
so  far  confined  to  two  dimensions  and  simple  substructures  with  benign  aspect  ratios, 
show  great  promise.  In  addition  to  Poisson's  equation,  membrane  and  shell  models 
have  been  tested.  So  far,  the  condition  numbers  have  always  been  less  than  4. 

The  following  theoretical  result  has  been  established  for  both  two  and  three  di- 
mensions. It  also  holds  for  elliptic  systems  such  as  those  of  linear  elasticity;  cf.  Smith 

[48]. 

Theorem  5.  The  operator  P  of  the  additive  algorithm  developed  by  Smith  satisfies 
the  estimate  k{P)  <  const.  The  constant  in  the  estimate  is  independent  of  h  as  well 
as  H. 
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7.  Optimal  Iterative  Refinement  Methods.  In  our  study  of  iterative  refine- 
ment methods,  we  consider  problems  on  a  special  kind  of  composite  finite  element 
triangulations.  We  begin  by  introducing  a  relatively  coarse  triangulation  of  Q^''  =  J7, 
and  denote  the  corresponding  space  of  finite  element  functions  by  V^^ .  We  can  think 
of  this  space  as  having  a  uniform  (or  relatively  uniform)  mesh  size  hi.  Let  fi*^^ 
be  a  subregion  where  we  wish  to  increase  the  resolution.  We  do  so  by  subdividing 
the  elements  and  introducing  an  additional  finite  element  space  V'^'^.  This  space  is 
constructed  quite  similarly  to  the  previous  one  and  it  contains  V^^  H-^ol^z)  as  a 
subspace.  We  assure  that  the  resulting  composite  space  V'^^  +  V'^^  is  conforming  by 
having  the  functions  of  V^'^  vanish  on  d^^^\  We  repeat  this  process  by  selecting  a 
subregion  i}^^^  of  17^^'  and  introducing  a  further  refinement  of  the  mesh  and  the  finite 
element  space,  etc.  We  denote  the  resulting  nested  subregions  and  subspaces  by  Q'*' 
and  V^' ,  respectively.  Throughout,  we  assume  that 

and  that 

V'''-'  nH^{n^'^)C  V'"'  C  ^o(^^'^)>      i  =  2,...,k. 
The  composite  finite  element  space,  on  the  repeatedly  refined  mesh,  is 

V^  =  f'''  +  v'^^  +  . . .  +  V'^'' 

The  finite  element  models  on  composite  meshes  is  thus  systematically  constructed 
by  introducing  a  basic  finite  element  approximation  on  the  entire  region  and  then 
selecting  subregions,  and  subregions  of  subregions  etc.,  where  the  finite  element  model 
is  further  refined  in  order  to  gain  higher  accuracy.  In  each  iteration  of  the  iterative 
refinement  methods,  problems  representing  finite  element  models  on  the  original  region 
and  the  subregions,  prior  to  further  refinement,  are  solved.  As  always,  an  additive 
algorithm  is  defined  by  specifying  the  subspaces  V,,  or  alternatively,  the  projections 

In  order  to  prove  the  results  given  at  the  end  of  this  section,  we  need  some  ad- 
ditional technical  assumptions.  To  make  our  proofs  work,  we  cannot  allow  the  sets 
fi,_i  \fi,  to  become  arbitrarily  thin  in  comparison  with  the  diameter  of  n,_i.  We  also 
assume  that  the  area  of  any  triangle  on  level  i  can  be  bounded  by  const.  q'~^  times 
the  area  of  the  triangle  on  level  j  of  which  it  is  a  part.  Here  g  is  a  constant  <  1  and 
j  <  I. 

The  fundamental  building  blocks  of  our  algorithms  are  the  projections  P/,  j  = 
i  —  l.z,  onto  the  spaces  V^'*^  n.ff^o(^t)-  ^^  ^^^^  ^^^^  iO  =  z  —  1,  we  solve  a  problem  on 
f2,  with  a  coarser  mesh  than  if  V'^'  were  used.  The  projection  P^  is  defined  in  terms 
of  the  unique  element  of  V^^  n  Hq{Q,),  which  satisfies 

(28)  aiP^VH,4>h)  =  aivH,<i>h),  'i4>h  G  V^' r\  Hli^i). 

We  consider  two  different  algorithms  and  distinguish  between  them  by  using  su- 
perscripts. The  perhaps  most  natural  algorithm  uses  the  projections  P^'  =  Pj.  The 
condition  number  of  this  algorithm  can  grow  as  fast  as  linearly  with  k.  By  the  stan- 
dard argument,  it  is  easy  to  show  that  the  eigenvalues  of  P^^)  are  bounded  from  above 
by  k.  This  bound  is  attained  if  V'"  n  H^{Qk)  is  not  empty,  i.e.  when  the  mesh  size  hi 
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is  fine  enough.  Any  function  in  this  space  belongs  to  V^' ,  i  =  1,2,...,/:,  and  is  exactly 
reproduced  by  each  of  the  projection  operators.  It  is  therefore  an  eigenfunction  of  P^^' 
with  the  eigenvalue  k.  Similarly,  any  function  which  belongs  to  V^^  f]  Hq{^i  \  O2) 
is  an  eigenfunction  with  eigenvalue  1.  We  have  shown  in  [56],  that  the  eigenvalues  of 
pO  are  bounded  from  below  by  a  constant.  The  condition  number  of  P^^)  is  therefore 
of  order  k. 

We  are  principally  interested  in  the  additive  algorithm  defined  by  the  projections 
pj  '  =  P*  —  P^'^j  ,  i  <  k  —  1,  and  P^  —  P^.  It  easy  to  show  that  these  operators  are 
projections  and  that  the  composite  finite  element  space  V^  is  the  direct  sum  of  the 
corresponding  subspaces  V,-  '.  The  following  results  have  recently  been  demonstrated 
in  Dryja  and  Widlund  [27]  and  Widlund  [56].  The  latter  paper  also  contains  a  similar 
result  for  a  multiplicative  algorithm. 

Theorem  6.  The  condition  number  of  P^^^  is  uniformly  bounded  by  a  constant. 
The  condition  number  of  P^^^  grows  at  most  linearly  with  the  number  of  refinement 
levels. 

We  note,  that  Tarek  Mathew  [40]  has  developed  iterative  refinement  methods 
for  Raviart-Thomas  finite  element  methods  and  obtained  bounds  for  their  rate  of 
convergence. 
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